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ABSTRACT 


A  least  squared  error  algorithm  is  presented  for  fitting  a  piecewise  linear  function  to 
observed  data,  when  the  slope  of  the  function  is  allowed  to  change  only  at  specified  points. 
The  algorithm  can  be  used  to  estimate  piecewise  linear  trends  in  data,  where  the  locations  of 
possible  changes  in  trend  are  known.  Furthermore,  it  cam  be  used  to  reduce  large  quantities  of 
data  to  manageable  sizes,  since  J  coordinates  are  input  and  L  4C  J  coordinates  are  output.  In 
addition,  the  algorithm  is  robust  to  problems  of  data  dropout  provided  it  is  used  cautiously. 
An  example  of  a  possible  application  is  the  fitting  of  a  linear  segment  bathythermal  profile 
(temperature  vs.  depth)  to  a  large  quantity  of  data. 


RESUME 

On  presente  un  algorithme  de  correction  d’erreur  par  la  methode  des  moindres  carres  qui 
permet  d’ajuster  une  fonction  lineaire  par  morceaux  a  un  ensemble  de  donnees  d’observation, 
quand  la  pente  de  la  fonction  ne  peut  varier  qu’a  des  points  precis.  Cet  algorithme  peut  servir 
4  estimer  les  tendances  lineaires  par  morceaux  des  donnees,  a  la  condition  que  soient  con- 
nus  les  points  des  variations  possibles  des  tendances.  Q  peut  aussi  etre  employe  pour  reduire 
de  grands  ensembles  de  donnees  a  des  ensembles  de  tailles  abordables,  puisque  pour  J  coor- 
donnees  d’entree  il  produit  L  J  coordonnees  de  sortie.  Enfin,  l’algorithme  est  insensible  aux 
problemes  des  pertes  de  donnees,  pourvu  qu’on  s’en  serve  avec  precaution.  Une  de  ses  applica¬ 
tions  possibles  est  l’ajustement  d’un  profil  bathythermique  lineaire  par  segments  (temperature 
en  fonction  de  la  profondeur)  a  un  grand  ensemble  de  donnees. 
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1  Introduction 


A  least  squared  error  algorithm  is  presented  for  fitting  a  piecewise  linear  function  to 
observed  data.  The  slope  of  the  function  is  allowed  to  change  only  at  specified  points.  This 
algorithm  can  be  used  to  estimate  piecewise  linear  trends  in  data,  where  the  locations  of  possible 
changes  in  trend  are  known.  A  possible  application  in  underwater  acoustics  is  the  fitting  of  a 
linear  segment  bathythermal  profile  (temperature  vs.  depth)  to  a  large  quantity  of  observed 
data. 


1.1  Statement  of  the  Problem 

The  observed  data  consist  of  J  coordinate  pairs  (xy, yy } , J  =  1  ,...J  where  xy  is  the 
independent  variable  and  y y  is  the  dependent  variable.  The  objective  is  to  fit  a  piecewise  linear 
function  to  this  data.  The  function  is  allowed  to  change  slope  at  x  =  uj,  l  =  1, . .  .L.  The 
values  are  given,  in  monotonic  increasing  order  <  uj  <  . . .  <  u i.  The  fitted  curve  can 
be  described  by  its  coordinates  (uj,  vj},/  =  1  ,...L  where  L  is  much  smaller  than  J.  The  fitted 
curve  should  be  one  which  minimizes  the  sum  of  the  squared  errors  between  the  observed  data 
and  the  fitted  curve,  with  allowance  for  possible  weighting  of  the  errors.  The  piecewise  linear 
function  /(x)  takes  the  following  form: 

X  ““  Uj 

/(x)  =  VI  + - [vi+1  -  Vi]  ;  U|  <  x  <  U|+1  (1.1) 

UJ+1  -  uj 


1.2  The  Least  Squares  Solution 

The  domain  of  x  is  divided  into  the  following  L  fixed  intervals  /j: 

h  =  (-oo.uj)  U  [U!  ,  u2) 

I2  =  [u2  ,  Uj) 

h  =  [«j  ,  u«+i) 

h  =  [«l  .  °°) 

An  additional  point  u^+i  is  defined  as  follows: 

1  r  i 

«Z,+  1  =  «£  +  ~  “‘I  ' 


(1.2) 


(1.3) 


1 


'A 


The  function  /(z)  is  rewritten  as  a  function  of  a  known  constant  /o  plus  an  unknown 
intercept  ao  plus  the  unknown  function  increments  aj  =  vi+j  -  vj: 

/(X)  =  /o  +  «o+  E  +  «<(.)  i  1  £  %  •  C-4) 

The  error  cy  associated  with  the  data  point  {zy,  yy}  is  defined  as  the  difference  between  the 
fitted  function  /(zy)  and  the  observed  value  yy: 


/(zy)  +  ey=  yy  ;  j  =  l,...J 

l(i)- l  /  _  \ 

«o+  ga*  +«w(ui(^1!'(;|(.|)+<i=  W*  i 


(1.5) 


This  is  a  system  of  J  linear  equations  in  L  +  1  unknowns  a/,  which  may  be  written  in  matrix 
notation  as  follows: 


Bjx(£+l)a(L+l)xl  +  e/xi  =yjxi  -  /o'ljxi 
where  row  j  of  the  matrix  B  is  as  follows: 


(1.6) 


1(0),  l(i),  ...  l(i(y)-i),  (  %i  ■  U<W  ),  0,...,0  (1.7) 

L  V  “»(;)+!  "  u«U)/ 


and  where  the  vector  of  unknowns  a  is  as  follows: 

aT  =  [ao,al,...,0£(] 


(1.8) 


The  least  squares  solution  a  minimizes  ere  =  £y=i  ey .  The  solution  is  obtained  by  setting 
the  gradient  of  this  function  with  respect  to  a  equal  to  zero: 


V,  (ere)  =  0  . 

This  results  in  the  following  system  of  linear  equations: 


(1.9) 


[BTB](i+i)x(L+i)a(£+1)x1  =  [Br(y  -  /o  •  l)](L+i)xi  •  (110) 

The  inverse  of  the  matrix  BTB  exists  if  and  only  if  BrB  has  full  rank  L  +  1.  This 
condition  is  satisfied  if  the  zy  are  unique,  J  >  L  +  1  and  if  each  interval  It  contains  at  least  one 
interior  data  point  Zy.  If  empty  intervals  occur,  the  problem  can  be  decomposed  into  smaller 
subproblems  where  each  subproblem  contains  no  empty  intervals.  In  this  case  /(z)  on  an  empty 
interval  is  defined  to  be  the  straight  line  segment  joining  the  end  of  the  last  nonempty  interval 
with  the  beginning  of  the  next  nonempty  interval. 
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A  weighted  least  squares  solution  is  one  which  minimizes  J3/=i  wjej  w^ere  the  wj  are 
positive  weights.  The  weighted  solution  is  obtained  by  multiplying  row  j  of  matrix  B  and  the 
right  hand  side  yj  -  fo  by  the  constant  y/vij,  j  =  1 .J. 

For  numerical  stability  the  constant  /q  is  chosen  to  be  the  mean  of  the  observed  values 
yj .  The  condition  number  of  BrB  is  best  if  the  data  points  are  located  new  the  centers  of  the 
intervals. 
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2  The  Least  Squares  Algorithm 

2.1  Initialization 

The  algorithm  begins  with  a  first  pass  at  the  observed  data  {xj,  y;)  to  compute  the  mean: 

/o  =  y  =  7  II V;  (2-1) 

J  i«i 

and  to  count  the  number  of  data  points  it]  in  each  interval  7j,  l  = 

2.2  Weighted  Errors 

The  error  weighting  function  may  be  chosen  as  required,  depending  on  the  application. 
For  illustration,  it  will  be  assumed  that  all  data  points  within  a  given  interval  should  be  equally 


weighted.  It  will  also  be  assumed  that  the  average  error  for  each  nonempty  interval  should  be 
equally  weighted.  The  following  weighting  function  reflects  these  assumptions: 


1 


c  + 


z;  G  h{i)  > 


(2.2) 


where  c  >  0  is  a  constant  which  may  be  used  to  modify  the  effect  of  segments  which  contain 
very  few  points.  In  the  examples  which  follow,  c  equals  10. 


2.3  Subproblem  1 


The  first  subproblem  is  numbered  n  =  1.  The  values  rq  are  scanned  to  find  the  first  set  of 
nonempty  adjoining  intervals  fix  =  Subproblem  n  has  L(n)  intervals 

with  L(n)  +  1  unknowns. 


2.3.1  Calculation  of  the  Least  Squares  Matrix 


The  coefficients  of  the  least  squares  system  in  Equation  1.10  are  obtained  by  processing 
the  data  set  {xj,y}}  point  by  point  without  bringing  the  entire  data  set  into  the  computer 
memory.  It  is  assumed  that  the  data  are  in  random  order;  therefore  the  entire  data  set  must  be 
scanned  in  order  to  solve  a  subproblem.  The  matrix  C  =  BrB  and  the  right  hand  side  vector 
d  =  BT[y  -  /o  •  1]  are  initially  set  to  zero: 


Cim=  0  ;  i  =  l,...L(n)  +  1  ;  m  =  1, . .  .L(n)  +  1 
di=  0  ;  i  =  1,.  ..L(n)  +  1  . 


(2.3) 


The  data  point  {xj,  y;}  is  contained  in  interval  I^y  If  fjy)  <f.  Qn  the  data  point  is  ignored. 
If  C  fin  then  the  following  extended  row  vector  b  is  defined  with  dimension  L(n)  +  2: 


b  = 


1(o) 


,  id),  itomw).  •  0 . °-  •  (2-4> 


The  jth  data  point  then  makes  the  following  contribution  to  C  and  d: 

C\m-  Cim  +  wjbibm  ;  »  =  l,...L(n)  +  l  ;  m  =  1, . .  .L[n)  +  1 

dj  =  di  +  «>/Mi,(n)+ 2  !  *  =  1,  •  •  •  L\n)  +  1 
where  the  prime  denotes  an  update  to  the  coefficient  values. 


(2.5) 


2.3.2  Solution  of  the  Least  Squares  System 

The  system  of  linear  equations  is 

Ca  =  d  (2.6) 

Since  the  matrix  C  is  symmetric,  an  efficient  inversion  algorithm  could  be  invoked.  However 
there  are  cases  where  C  may  be  ill-conditioned  if  the  data  are  sparse  and  if  some  of  the  data 
points  {xj,  y;)  are  very  close  together,  A  modified  gaussian  elimination  technique  has  been 
selected  instead.  This  technique  is  based  on  the  Simplex  Method  of  linear  programming  [l] .  If 
the  matrix  C  has  rank  M  <  L(n)  +  1,  then  the  algorithm  assigns  arbitrary  values  to  L(n)  +  l- M 
variables.  The  preferred  arbitrary  variables  are  ajvr, . . . ,  aL{n)-  Arbitrary  default  values  for  all 
variables  should  be  chosen  to  suit  the  application.  For  bathythermal  curve  fitting  the  default 
values  may  be  set  to  zero,  or  they  may  be  chosen  based  on  historical  data. 


2.3.3  Assignment  of  Coordinates  of  f(x) 

After  solution  of  the  first  subproblem  n  =  1  the  solution  of  /(x)  is  determined  for  x  < 
u*(i)+£(i).  The  coordinates  vj  are  calculated  as  follows: 


( 


vi  =  /(«/)  = 


fo  +  ao  +  a\ 


f  ui  -  Uk(i)  \ 
\vkW+i-uk{l)J 


i-i 

/o  +  ao+  E  al+m-i(l) 

m=i(l) 


l  <  A:(l) 


Jfc(l)  <  /  <  Jt(l)  +  £(1)  . 


(2.7) 


This  assumes  that  if  the  first  intervals  contain  no  data  points,  the  initial  coordinates  are  ob¬ 
tained  by  linear  extrapolation  of  the  first  nonempty  segment;  these  coordinates  could  also  be 
specified  arbitrarily,  depending  on  the  application. 


2.4  Termination  Criteria 

At  the  end  of  subproblem  n,  one  of  the  following  cases  will  apply: 

Case  I.  If  k(n )  +  L(n)  =  L,  then  /(x)  has  been  determined  for  all  x,  and  there  are  no 
more  subproblems  to  solve. 

Case  II.  If  fc(n)  +  L(n)  <  L,  and  the  remaining  intervals  are  empty,  then  /(x)  may  be 
extrapolated  linearly  from  the  last  segment: 

vi  =  /(uj) 

- Ul  -  - )  ;  /  >  k[n)  +  L(n) 

Uk(n)+L(n)  Uk(n)+L{n)~  1  J 

(2-8) 

or  the  coordinates  may  be  specified  arbitrarily.  The  function  /(x)  is  then  completely  defined 
and  there  are  no  more  subproblems  to  solve. 

Case  HI.  If  k{n)  +  L(n)  <  L  and  there  are  non-empty  intervals  remaining,  i.e. 

E  ,  (2.9) 

i=*(n)+L(r»)+l 

then  at  least  one  more  subproblem  must  be  solved. 


fc(n)+L(n)-l 

=  fo  +  OO  +  E  al+m-*(n)  +  <*£,{„) 

m=k(n) 
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2.5  Subproblem  n 


The  nth  subproblem  begins  with  the  identification  of  the  next  set  of  nonempty  adjoining 
intervals  f2„  =  {h(n)>- ■ -h{n)-i+L(n)}-  The  new  least  squares  matrix  C  and  vector  d  are 
computed  and  the  system  is  solved  for  the  new  unknowns  a.  At  the  end  of  subproblem  n  (n  >  1) 
the  solution  of  /  (i)  is  determined  for 

U/k(n)  <X<  Uk(n)+L(n)  .  (2-10) 


The  coordinates  t/j  are  calculated  as  follows: 
(  fo  +  ao 


vi  =  /(uj)  =  < 


/-l 

fo  +  <k)  +  £  al+m-k(n) 

m=k(n) 


!  l  =  k(n) 

;  k(n)  <  l  <  k[n)  +  L[n)  . 


(2.11) 


The  intermediate  coordinates  between  problem  n  -  1  and  problem  n  are  defined  as  follows: 


vi  =  V*(n-l)+£(n-l)  + 


vk(n)  Vib(n-1)+X.(n-1) 
uk(r»)  ~  Uk(n-1)+L(n-1) 


(«!-«k(n-l)+L(n-l)) 


A:(n-1)+L(n-1)  <  /  <  k(n) 
(2.12) 


The  algorithm  then  determines  which  of  the  three  termination  cases  applies  at  the  end 
of  subproblem  n.  If  Case  I  or  Case  II  applies,  the  algorithm  terminates.  If  Case  III  applies, 
another  subproblem  must  be  solved. 


3  Examples 


A  few  examples  have  been  chosen  to  demonstrate  the  curve  fitting  procedure  in  extreme 
cases.  The  bathythermal  data  are  from  shallow  water.  Appendix  A  contains  a  listing  of 
FORTRAN  program  FITBT  which  specifies  the  required  values  of  ui  and  default  values  for  the 
variables  aj.  Appendix  B  contains  a  listing  of  FORTRAN  subroutine  LSQPL  which  implements 
the  least-squares  algorithm. 


3.1  Example  1:  A  Well-Behaved  Problem 

Figure  3.1  shows  a  well-behaved  example  in  which  there  are  many  data  points.  The  fixed 
depths  U(  are  indicated  by  dashed  lines.  Since  all  of  the  depth  intervals  contain  data  points, 
there  is  only  one  problem  to  solve.  Figure  3.2  shows  the  fitted  piecewise  linear  curve. 


3.2  Example  2:  Problem  Decomposes  into  Two  Subproblems 

Figure  3.3  shows  an  example  in  which  there  are  several  empty  intervals.  The  problem 
decomposes  into  two  independent  subproblems,  as  shown  in  Figure  3.4.  The  intermediate  co¬ 
ordinates  are  defined  by  a  linear  segment  joining  the  two  independent  solutions.  The  first 
coordinates  are  obtained  by  linear  extrapolation  of  the  first  nonempty  interval.  The  last  coor¬ 
dinates  are  obtained  by  linear  extrapolation  of  the  last  nonempty  interval. 


3.3  Example  3:  An  Example  with  Sparse  Data 

Figure  3.5  is  an  extreme  example  in  which  the  number  of  unknowns  exceeds  the  number 
of  data  points.  One  solution,  shown  in  Figure  3.6,  sets  the  free  variables  to  zero;  that  is,  the 
temperature  increments  in  the  associated  intervals  (4  and  10)  are  zero.  A  second  solution, 
shown  in  Figure  3.7,  chooses  the  free  variable  values  based  on  historical  temperature  gradient 
data,  Figure  3.8.  Appendix  C  contains  the  corresponding  input  and  output  files  for  example  3. 
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Figure  3.1:  Example  1.  A  well-behaved  example  with  more  data  points  than  unknowns.  The 
fixed  depths  tij  are  indicated  by  dashed  lines. 


Observed  Data 


Temperature  y  (#C) 

Figure  3.5:  Example  3.  Sparse  Data 


q-q  Curve  fit 


Temperature  y  (°C) 

Figure  3.6:  Example  3  -  Solution  a.  Solution  of  a  degenerate  problem  allows  free  function 
increment  variables  to  be  set  to  arbitrary  user-defined  values.  In  this  case  the  free  variables 
(temperature  increments  in  intervals  4  and  10)  are  defined  to  be  zero. 


4  Summary 

An  algorithm  has  been  presented  for  obtaining  a  weighted  least-squares  fit  of  a  piecewise 
linear  model  to  observed  data,  when  the  linear  segment  domains  are  specified.  This  algorithm 
may  be  used  to  reduce  large  quantities  of  data  to  more  manageable  size,  since  J  coordinates 
are  input  and  L  «  J  coordinates  are  output.  The  algorithm  is  reasonably  robust  to  problems 
of  data  dropout,  but  should  be  used  cautiously.  Caution  is  recommended  because  sparse  data 
in  one  segment  may  have  undue  influence  on  the  curve  fit  unless  suitable  weighting  functions 
are  defined. 


Appendix  A  Listing  of  FORTRAN 
Program  FITBT 


PROGRAM  FITBT 

C  Fit  a  Piecewise  Linear  Curve  to  Bathythermal  Data  (Depth  vs.  Temperature) 
C  or  to  Sound  Profile  Data  (Depth  vs.  Sound  Speed) 

C  using  Least  Squares  Algorithm  LSQPL  and  Linear  Equation  solver  SOLV. 

C 

C  Up  to  100  unknowns  (99  Linear  Segments) 

C*****************^***:|^:|^**:|^**************!|t************:|t******************** 

C  Input  observed  data  from  disk  file  LSQIN.DAT 
C  J  (Number  of  data  points) 

C  Yl,  XI  (TEMPERATURE  OR  SPEED  Tl.  DEPTH  XI) 

C  . 

C  YJ.  XJ  (TEMPERATURE  OR  SPEED  YJ.  DEPTH  XJ) 

C  No  internal  memory  limit  on  the  number  of  data  points. 

C*****4^**************4^i|^***:|^*********************************************** 

C  Optional  input  of  historical  data  (up  to  100  points)  from  disk  file 
C  SPRF0.DAT 

C  used  to  define  the  number  of  segments,  L,  the  fixed  depths  U(I), 

C  and  default  values  for  function  increment  variables  DELV(I) 

C  (22-CHARACTER  IDENTIFIER) ,  L 

C  U(1).TEMP(1),SPEED(1)  (HISTORICAL  DEPTH,  TEMPERATURE.  SOUND  SPEED) 

C  . 

C  U(L) , TEMP (L) , SPEED (L) 

C***********************^***************************************************** 

C  Output  fitted  data  on  disk  file  FITBT.DAT 
C  L  (Number  of  fitted  points  L«J) 

C  V(l),  U(l)  (TEMPERATURE  OR  SPEED  V(l),  DEPTH  U(l)) 

C  . 

C  V(L) ,  U(L)  (TEMPERATURE  OR  SPEED  V(L) ,  DEPTH  U(L)) 

C***************** ************ ************************************************* 

CHARACTER  UNITS,  ID(22) 

DIMENSION  UU(101) 

DIMENSION  U(101) ,V(101) .DELV(IOO) 

DIMENSION  A(101) ,B(102) ,C(101,102) ,NL(100) . INDEX(lOl) , JNDEX(lOl) 
DIMENSION  TEMP (101) .SPEED (101) 

C  STANDARD  DEPTHS  IN  METERS  FOR  DEEP  WATER  BATHYTHERMAL  PROFILE 


C  TWO  ADDITIONAL  DEPTHS  (62.5,87.5)  HAVE  BEEN  ADDED  FOR  SHALLOW  WATER 
DATA  DU/ 

1  0.0,10. ,20. ,30. .50. ,62.5,75. ,87.5,100. .125. ,150. .200. ,250. .300. . 

2  400. .500. ,600. .700. ,800. ,900. ,1000. ,1100. ,1200. .1300. ,1400. . 

3  1500 . . 1750 . , 2000 . , 2500 . , 3000 . , 4000 . , 5000 . , 69*0 . / 

C  CONST  -  CONSTANT  TO  BE  USED  IN  WEIGHTED  ERROR  FUNCTION.  CONST>»0. 

C  THE  WEIGHT  ON  THE  ERROR  ASSOCIATED  WITH  OBSERVATION  JJ 

C  -  1/ (CONST  +  NO.  OF  POINTS  IN  THE  SAME  INTERVAL  AS  DATA  POINT  JJ) 

CONST  •  10. 

LMAX  -  100 
LMAX1-  LMAX+1 
LMAX2-  LMAX +2 
L  -  30 

DO  400  I-1.LMAX1 
400  U(I)«UU(I) 

WRITE (5, 500) 

500  FORMAT ( IX’  PROGRAM  FITBT7 

1  IX ‘VERSION  850708  -  DREA  -  B . A . TRENHOLM 7 

2  IX ‘OBSERVED  DATA  FROM  DISK  FILE  LSQIN.DAT:  7 

3  IX*  IN  ENGLISH  UNITS (0)  OR  METRIC (1)  ?  :  ’,$) 

READ (5 , *) IFTM 

UNITS-  ‘ft* 

IF(IFTM . NE . 0)UNITS-  ’m’ 

WRITE (5, 501) 

601  FORMAT ( 

1  1X‘  TEMPERATURE  DATA(O)  OR  SOUND  SPEED  DATA(l)  ?  :  ’,$) 

READ(5,*)IBTVEL 

C  ***  DEFINE  VALUES  DELV(I)  TO  CONTROL  DEFAULT  VALUES  FOR 
C  FUNCTION  INCREASE  IN  LAYER  (I)  IF  EQUATIONS  ARE  UNDERDETERMINED. 

U (L+l ) -U (L) + (U (L) -U( 1 ) ) /FLOAT (L- 1 ) 

DO  502  I-l.L 

C  DEFAULT  TO  ZERO  TEMPERATURE  GRADIENT 

DELV(I)-0 .0 

C  DEFAULT  TO  SOUND  SPEED  PRESSURE  GRADIENT 

IF(IBTVEL.NE.0)DELV(I)«0.017*(U(I+1)-U(I)) 

502  CONTINUE 


INFTM-1 

C  ***  OPTION  TO  INPUT  HISTORICAL  DATA 
WRITE(5,504) 

504  FORMAT  ( 

1  lX'READ  HISTORICAL  DATA  FROM  FILE  SPRFO.DAT  (O-NO,  1-YES) :  \$) 
READ(5 ,*)IHIST 
IF ( IHIST . EQ . 0) GO  TO  550 
WRITE(5,505) 

605  FORMAT ( IX *  ENGLISH  UNITS  (0)  OR  METRIC  (1)  ?  :  \$) 

READ (5 ,*) INFTM 

OPEN (UNIT-23, DEVICE- ’DSK’ .FILE-’ SPRFO.DAT’ , ACCESS- ’SEqiN’ , 

1  MODE-* ASCII’) 

READ (23 .2300) ID ,L 
2300  FORMAT (22A1 , 15) 

IF (L . GT . LMAX) WRITE (5 , 506) L , LMAX 

506  F0RMAT(1X'L-’I4, ’ .  MORE  THAN  *  14 .’  SEGMENTS !’ ) 

IF (L.GT. LMAX) STOP 

DO  507  I-l.L 

READ(23,*)U(I) ,TEMP(I) ,SPEED(I) 

IF(I.Eq.l)GO  TO  507 
IF (U(I).LE.U(I-1))WRITE(5. 608)1, U(I) 

508  FORMAT ( IX ’ DEPTH ’ 14 , ’-’Ell .6, '  IS  NOT  IN  INCREASING  ORDER’) 

IF(U(I).LE.U(I-1))ST0P 
DELV(I-1)-TEMP(I)-TEMP(I-1) 
IF(IBTVEL.NE.0)DELV(I-1)-SPEED(I)-SPEED(I-1) 

507  CONTINUE 

U(L+1)-U(L)+(U(L) -U(l) )/FLOAT(L-l) 
DELV(L)-DELV(L-1)*(U(L+1)-U(L))/(U(L)-U(L-1)) 


C  ***  CHECK  TO  SEE  IF  UNIT  CONVERSION  IS  REQUIRED 
650  IF ( INFTM . EQ . IFTM) GO  TO  559 

IF ( INFTM . EQ . 0) GO  TO  540 

C  ***  CONVERT  DEPTHS  AND  DEFAULT  VARIABLES  TO  ENGLISH  UNITS 
DO  639  1-1 ,L 
U(I)-U(I)*3. 28084 

IF ( IBTVEL . NE . 0) DELV ( I ) -DELV ( I ) *  3 . 28084 
IF(IBTVEL. EQ .0)DELV(I)-DELV(I) *9 . /5 . 

639  CONTINUE 
GO  TO  559 


C  *«■*  CONVERT  DEPTHS  AND  DEFAULT  VARIABLES  TO  METRIC  UNITS 


540  DO  540  I-l.L 

U(I)-U(I)/3. 28084 

IF (IBTVEL . NE . 0)DELV (I) -DELV(I) /3 . 28084 
IF ( IBTVEL . EQ . 0 ) DELV ( I ) -DELV (I ) * 5 . /9 . 

549  CONTINUE 

C  ***  DEFINE  MAXIMUM  DEPTH  IN  FITTED  PROFILE 
559  INFTM-IFTM 

WRITE (5, 503) UNITS 

603  F0RMAT(1X* INPUT  WATER  DEPTH  (’A2, ’)  :  *,$) 

READ(5 ,*)WD 

DO  600  1*1 ,L 
IF(WD.LE.U(I))GO  TO  601 

600  CONTINUE 
I  -  L+l 

601  L  -  I 
U(L)  *  WD 

C  ***  IF  THERE  ARE  NO  DATA  POINTS  IN  THE  FIRST  SEGMENTS,  USE  LINEAR 
C  EXTRAPOLATION  TO  DEFINE  COORDINATE  VALUES 

ISETA-0 

C  ***  IF  THERE  ARE  NO  DATA  POINTS  IN  THE  LAST  SEGMENTS,  USE  LINEAR 

C  EXTRAPOLATION  TO  DEFINE  COORDINATE  VALUES 

ISETB-0 

C  *  HOWEVER,  IF  HISTORICAL  DATA  ARE  AVAILABLE,  USE  DELV(I)  TO 
C  OBTAIN  COORDINATE  VALUES 

IF ( IHIST . NE . 0) ISETB* 1 

C  ***  IF  THE  CURVE  FIT  PROBLEM  IS  UNDERDEFINED  (MORE  VARIABLES  THAN 
C  DATA  POINTS)  USE  DELV(I)  TO  DEFINE  FREE  VARIABLES. 

ISET-1 

C  ***  OPTIONAL  DEBUG  OUTPUT  TO  TERMINAL 
WRITE (5, 5029) 

5029  FORMAT (IX ’OPTIONAL  DEBUG  OUTPUT  TO  TERMINAL  (0*N0,  1-YES)  ?  * 
READ(5 , *) IBUG 

C  ***  CALL  SUBROUTINE  TO  SET  UP  EQUATIONS  AND  SOLVE  FOR  UNKNOWNS 
CALL  LSqPL (LMAX , LMAX1 , LMAX2 , L . A , B , C . NL . INDEX , JNDEX , ISET , 

1  I SETA , ISETB , U , V , DELV , N , IICK , CONST , IBUG) 


WRITE(5 ,5030) (I ,NL(I) , 1*1 ,L) 

5030  FORMAT (IX 1  LAYER  # ’ 13 . '  CONTAINS '  15 , '  POINTS ' ) 

WRITE (5 ,5010) N 
5010  FORMAT ( 

2  1X13, ’ SUBPROBLEMS  IDENTIFIED*/) 

IF(IICK.NE.0)WRITE(5 ,5020) IICK 
5020  FORMAT (1X13 'SUBPROBLEMS  WERE  INCONSISTENT  !') 

IF(IICK.EQ.O) 

1  OPEN (UNIT-22, DEVICE- 'DSK* .FILE- ’FITBT.DAT' , ACCESS- ' SEQOUT * , 

2  MODE- 'ASCII') 

IF ( I ICK . EQ . 0) WRITE ( 5 . 5040) 

5040  F0RMAT(1X' COORDINATES  OF  THE  FITTED  CURVE  ARE  ON  FILE  FITBT.DAT') 

IF(IICK.EQ.0)WRITE(22,*)L 
DO  220  I-l.L 

IF(IICK.EQ.0)WRITE(22,*)V(I),U(I) 

220  CONTINUE 


IF(IICK.EQ.O)  CLOSE (UNIT-22) 


Appendix  B  Listing  of  FORTRAN 
Subroutine  LSQPL 


SUBROUTINE  LSQPL (LMAX , LMAX1 , LMAX2 , L . A , B . C . NL . INDEX , JNDEX . ISET , 

1  ISETA . ISETB .U.V.DELV.N, IICK, CONST , IBUG) 

C 

C  B . A . TRENHOLM ,  D.R.E.A.,  JULY  1985. 

C  LEAST-SqUARES  FIT  OF  PIECEWISE  LINEAR  FUNCTION 

C************************** ******* ************************************* 

C  INPUT  VARIABLES: 

C  LMAX  »  MAXIMUM  NO.  OF  LINEAR  SEGMENTS 
C  LMAX1  -  LMAX+1 

C  LMAX2  -  LMAX +2 

C  L  -  NO.  OF  LINEAR  SEGMENTS  TO  BE  FITTED  TO  DATA  (L.LE.LMAX) 

C  U(K)  -  FIXED  X  COORDINATES  U(1)<U(2)<. . .<U(L)  DEFINING  SEGMENTS 

C  A.B.C.NL.INDEX. JNDEX  ARE  WORKING  ARRAYS 
C  ISET  -  0  IF  FREE  VARIABLES  ARE  SET  TO  ZERO 

C  »  1  IF  FREE  VARIABLES  DEPEND  ON  THE  INPUT  VALUES  DELV(K) 

C  I.E.  FREE  VARIABLE  A(H-K)-DELV(KMIN+K-1) 

C  (FREE  VARIABLES  OCCUR  WITHIN  A  SUBPROBLEM  IF  THERE  ARE 

C  MORE  UNKNOWNS  THAN  DATA  POINTS.) 

C  ISETA  -  0  LINEAR  EXTRAPOLATION  WILL  BE  USED  TO  DEFINE  FIRST 
C  COORDINATES  IF  THE  FIRST  SEGMENTS  CONTAIN  NO  DATA  POINTS. 

C  -  1  DELV(I)  WILL  BE  USED  TO  DEFINE  FIRST  COORDINATES 

C  IF  THE  FIRST  SEGMENTS  CONTAIN  NO  DATA  POINTS. 

C  ISETB  -  0  LINEAR  EXTRAPOLATION  WILL  BE  USED  TO  DEFINE  LAST 
C  COORDINATES  IF  THE  LAST  SEGMENTS  CONTAIN  NO  DATA  POINTS. 

C  -  1  DELV(I)  WILL  BE  USED  TO  DEFINE  LAST  COORDINATES 

C  IF  THE  LAST  SEGMENTS  CONTAIN  NO  DATA  POINTS. 

C  DELV(I)-  DEFAULT  VALUES  FOR  DIFFERENCES  V(I+1)-V(I),  I-1...L 
C  USE  OF  DELV  MAY  BE  RESTRICTED  BY  ISET . ISETA , ISETB . 

C  CONST  -  CONSTANT  TO  BE  USED  IN  THE  WEIGHTED  ERROR  FUNCTION.  CONST>*0. 

C  THE  WEIGHT  ASSIGNED  TO  THE  ERROR  ASSOCIATED  WITH  OBSERVATION  JJ 

C  -  1/ (CONST  +  NO.  OF  POINTS  IN  SAME  INTERVAL  AS  DATA  POINT  JJ) . 

C 

C  IBUG  -  0  TO  SUPPRESS  DEBUG  OUTPUT  TO  TERMINAL 
C  -  1  FOR  DEBUG  OUTPUT  TO  TERMINAL 

C 

C  INPUT  FROM  DISK  FILE  LSQIN.DAT: 


C  J  NUMBER  OF  DATA  POINTS 

C  Y1,X1  -  Y  AND  X  COORDINATES  OF  DATA  POINT  1 

C 

C  YJ.XJ  -  Y  AND  X  COORDINATES  OF  DATA  POINT  J 

C***^*^************************************************************** 

C  OUTPUT  VARIABLES: 

C  NLCK)  -  NUMBER  OF  DATA  POINTS  IN  INTERVAL  K,  WHERE 
C  INTERVAL  1  -  (-INFINITY, U(2)) 

C  INTERVAL  K  -  [U(K) ,U(K+1) ) 

C  INTERVAL  L  -  [U(L) ,+INFINITY) 

C  V(K)  -  VALUE  OF  FITTED  CURVE  AT  U(K) 

C  N  NUMBER  OF  SUBPROBLEMS  IDENTIFIED,  EACH  HAVING  CONSECUTIVE 

C  NONEMPTY  INTERVALS 

C  IICK  -  NUMBER  OF  INFEASIBLE  SUBPROBLEMS 

C**************************************** ***************************** 

C  INTERMEDIATE  WORKING  VARIABLES  FOR  NTH  SUBPROBLEM 

C  A(l)  -  INTERCEPT  AT  U(KMIN)  OFFSET  BY  FO,  I.E.  (V(KMIN)-FO  +A(1)) 

C  A(K)  -  FUNCTION  INCREMENT  IN  SEGMENT  (KMIN  +  K-2)  ,  K=*2, .  .  ,L(N)+1 

C  -  V(KMIN  +  K-l) -V(KMIN  +  K-2) 

C  B  »  EXTENDED  ROW  VECTOR  FOR  EQUATION  CORRESPONDING  TO  XJ.YJ 
C  C  AUGMENTED  LEAST  SQUARES  MATRIX  WITH  R.H.S.  VECTOR  IN  LAST  COL 

C  INDEX  -  INDEX  VECTOR  TO  MAP  SOLUTION  FROM  R.H.S.  OF  REDUCED  MATRIX 

C  JNDEX  »  INDEX  VECTOR  TO  IDENTIFY  FREE  VARIABLES  OF  DEGENERATE  MATRIX 

C  AFTER  MODIFIED  GAUSSIAN  ELIMINATION. 

C*************************************************************** ****** 

DIMENSION  A(LMAXl) ,B(LMAX2) ,C(LMAX1 ,LMAX2) , NL(LMAX) , INDEX (LMAX1) 
DIMENSION  JNDEX (LMAX1) 

DIMENSION  U(LMAXl) ,V(LMAX1) .DELV(LMAX) 

C  ADD  ONE  ADDITIONAL  POINT 

U(L+1)  -  U(L)  +(U(L)-U(1))/FL0AT(L-1) 

C  CLEAR  SOLUTION  VECTOR  AND  COUNTER 
DO  1  I-1,L 
NL(I)  «  0 

1  IF ( ISET . EQ.0)V(I)*0.0 


C  N-  NUMBER  OF  CURRENT  SUBPROBLEM 
N  -  0 
IICK  -  0 
KLAST  -  0 


KMINL  -  0 
KMAXL  -  0 


0PEN(UNIT=21, DEVICE- 'DSK' .FILE- 'LSQIN.DAT ’ .ACCESS- 'SEQIN' . 

1  MODE-’ ASCI I*) 

C  COMPUTE  AVERAGE  VALUE  OF  DATA  YJ.  COUNT  NUMBER  OF  POINTS  NL  IN  LAYER  LJ 
READ(21 , *) J 
IF (J.LE.O) RETURN 


FO  -  0.0 
DO  2100  I-l.J 
READ (21 , *)YJ ,XJ 
XI  -  I 

FO  -  (I-1)*(F0/XI)  +  YJ/XI 
LJ  -  1 

IF (X J . LT . U(2) ) GO  TO  2099 
LJ  -  2 

IF(L.EQ.2)G0  TO  2099 
DO  2098  K-3.L 
IF (XJ . LT . U (K) ) GO  TO  2099 
LJ  -  LJ  +1 

2098  CONTINUE 
LJ  -  L 

2099  NL(LJ)  -  NL(LJ)+1 

2100  CONTINUE 


2900  REWIND  21 
C  TEST  STATUS  OF  PROBLEM 

C  CASE  I  -  FUNCTION  IS  COMPLETELY  DEFINED 
IF (KLAST . GE . L) WRITE (5 , 51 1 )N 

511  FORMAT ( IX’ SUBPROBLEM ’ 13.  'TERMINATES  AS  CASE  I: ALL  VARS .DEFINED ' ) 

IF (KLAST . GE . L) CLOSE (UNIT-21 ) 

IF (KLAST. GE.L) RETURN 

C  FIND  NEXT  NONEMPTY  INTERVAL 

DO  100  I-KLAST+1 ,L 
KMIN  -  I 
KMAX  -  I 


21 
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IF (NL (KMIN) . EQ . 0) GO  TO  100 
C  NONEMPTY  INTERVAL  FOUND 

DO  98  K-KMIN.L 
IF(NL(K) .EQ.O)GQ  TO  104 
KMAX-K 

98  CONTINUE 

KMAX  -  L 
GO  TO  104 

100  CONTINUE 

C  CASE  II  -  REMAINING  INTERVALS  ARE  EMPTY 
IFCISETB .NE.O)GO  TO  102 

C  LINEAR  EXTRAPOLATION  OF  REMAINING  VARIABLES 
DO  101  I-KLAST+l.L 

101  V(I)-V(KMAXL)  +  DVDU*(U(I)-U(KMAXL)) 

WRITE(5,522)N 

522  FORMAT ( IX ’ SUBPROBLEM ’ 13 , '  TERMINATES  AS  CASE  Ha: EXTRAPOLATED’) 
CLOSE (UNIT-21) 

RETURN 

C  REMAINING  VARIABLES  CONTROLLED  BY  INPUT  VALUES  DELV(I) 

102  DO  103  I-KLAST+1 ,L 
V(I)-V(I-1)+DELV(I-1) 

103  CONTINUE 

WRITE(6,523)N 

523  FORMAT (IX’ SUBPROBLEM ’13,*  TERMINATES  AS  CASE  I lb: CONTROLLED*) 
CLOSE (UNIT-21) 

RETURN 

C  CASE  III  -  DEFINE  NEXT  SUBPROBLEM 

104  IF(N.GT.0)WRITE(5,533)N 

533  FORMAT (IX* SUBPROBLEM ’,13, ’  COMPLETE’) 

N  -  N  ♦  1 

C  ***  CLEAR  LEAST-SqUARES  MATRIX 
LN  -  KMAX-KMIN+1 
LN1  -  LN+1 
LN2  -  LN+2 

DO  105  I-l.LNl 


DO  105  K-1.LN2 
105  C(I,K)  -  0.0 

C  ***  READ  NUMBER  OF  POINTS  IN  DATA  FILE 
READ(21 ,*) J 

DO  260  JJ-l.J 
C  ***  READ  DATA  POINT  JJ 
READC21 ,*)YJ,XJ 

C  ***  CHECK  THAT  DATA  POINT  IS  CONTAINED  IN  INTERVALS  FOR  SUBPROBLEM  N 
IF (KMAX . LT . L . AND . (XJ . GE.U(KMIN) . AND.XJ.LT.U(KMAX+1)))G0  TO  250 
IF (KMAX . EQ . L . AND . (XJ.GE.U(KMIN)))GO  TO  250 
IF (KMIN . EQ . 1 . AND . X J . LT . U (2) ) 60  TO  250 
GO  TO  260 

C  ***  FIND  LAYER  LJ  CONTAINING  DATA  POINT 

250  LJ-1 

IF (X J . LT . U(2) ) GO  TO  258 
LJ-2 

IF(L.EQ.2)G0  TO  258 
DO  259  K-3.L 
IF (X J . LT . U (K) ) GO  TO  258 
LJ  -  LJ  +1 
259  CONTINUE 

LJ  -  L 

258  CONTINUE 

C  ***  EQUATION  FOR  DATA  POINT  XJ.YJ 
DO  251  I-1.LN2 

251  B(I)  -  0.0 
B(l)  -  1.0 

DO  252  I-l.LJ-KMIN 
BCI+1)  -  1.0 

252  CONTINUE 

BCLJ-KMIN+2)  -  (XJ-U(LJ))/(U(LJ+1)-U(LJ)) 

B(LN2)  ■  YJ-FO 

C  ***  WEIGHTED  ERROR  FUNCTION 

WJ  -  1 .0/ (CONST  ♦  FLOAT (NL(LJ))) 

C  ***  UPDATE  LEAST-SQUARES  MATRIX 
DO  263  I-l.LNl 


DO  253  K-1.LN2 

C(I ,K)  •  C(I ,K)  +  WJ*B(I)*B(K) 
253  CONTINUE 


260  CONTINUE 


C  ***  END  OF  DATA 


C  ***  DEBUG  PRINTOUT  OF  LEAST-SQUARES  MATRIX 
DO  5890  I-l.LNl 

IFCIBUG.LT. 0)WRITE(515891)I.(C(I1M).M-1.LN2) 
5991  FORMAT ( IX’ ROW 13.’:  ’100E12.5) 

5990  CONTINUE 


C  ***  DEFINE  DEFAULT  VALUES  FOR  FREE  VARIABLES 
IF(ISET.EQ.O)GQ  TO  262 
A(l)-0.0 
DO  261  I-2.LN1 

261  A(I)-DELV(KMIN+I-2) 

262  CONTINUE 


C  ***  SOLVE  LEAST-SqUARES  EQUATIONS 


CALL  SOLVCC . A . LMAX1 , LN1 , INDEX , JNDEX , ISET , ICK , IBUG) 


C  ***  COUNT  THE  NUMBER  OF  INFEASIBLE  PROBLEMS 
I ICK  -  I ICK  +  ICK 

IFCIBUG.NE.O. AND . ICK.EQ. 1) WRITE (5, 5992) N 
5992  FORMAT (IX ’SUBPROBLEM’ 13 , ’  IS  INCONSISTENT’) 


C  ***  DEFINE  COORDINATES  V(I)  OF  FITTED  FUNCTION  AT  U(I) 
IF(N.GT. 1)G0  TO  6002 


C  ***  FIRST  SUBPROBLEM  N-l,  DEFINE  COORDINATES  OF  LEADING  EMPTY  INTERVALS 
IFCKMIN.Eq. 1)G0  TO  6002 
IF(ISETA.NE.O)GO  TO  6000 
C  *  LINEAR  EXTRAPOLATION 
WRITE(5 ,513) 

513  FORMAT (IX’ INITIAL  COORDINATES  OBTAINED  FROM  EXTRAPOLATION’) 

DO  6001  I-l.KMIN-l 

V(I)  ■  FO  +  A ( 1 )  +A(2)*(U(I) -U(KMIN) )/(U(KMIN+l) -U(KMIN) ) 

6001  CONTINUE 


,*w  ■ 
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GO  TO  6002 

C  *  DETERMINED  BY  DELV(I) 

6000  WRITE(5,612) 

512  FORMAT (IX ’INITIAL  COORDINATES  CONTROLLED’) 

DO  6005  1-1 , KMIN- 1 
II-KMIN-I 

V(II)«V(II+1)-DELV(II) 

6006  CONTINUE 

C  ***  SUBPROBLEM  N,  COORDINATES  OBTAINED  FROM  LEAST-SQUARES  SOLUTION 

6002  DO  6004  I-KMIN.KMAX+1 
IF(I.GT.L)GO  TO  6004 
V(I)  -  FO  +  A(l) 

IF ( I . EQ . KMIN) GO  TO  6004 
DO  6004  K-KMIN.(I-l) 

V(I)  -  V(I)  +  AC2+K-KMIN) 

6004  CONTINUE 

C  ***  LINEAR  INTERPOLATION  BETWEEN  SUBPROBLEMS 
IF(N.EQ.1)GQ  TO  6009 
DVDU  «(V (KMIN) -VLAST)/(U (KMIN) -ULAST) 

DO  6007  I-KLAST+l.KMIN-1 
V(I)»VLAST+  DVDU* (U(I) -ULAST) 

6007  CONTINUE 

C  ***  STORE  ENDPOINT  OF  SUBPROBLEM  N 

6009  VLAST  -  V(KMAX+1) 

ULAST  -  U(KMAX+1) 

KLAST  -  KMAX+1 
KMINL  -  KMIN 
KMAXL  -  KMAX 

C  ***  PREPARE  FOR  LINEAR  EXTRAPOLATION  OF  REMAINING  COORDINATES 
IF (KLAST . EQ . L) GO  TO  2900 
DVDU  «  (VLAST-V(KMAX) )/ (ULAST-U(KMAX) ) 

GO  TO  2900 


Appendix  C  Listing  of  FORTRAN 
Subroutine  SOLV 


SUBROUTINE  SOLV (A , X.MAXN , N , INDX , JNDX , ISET , ICK , IBUG) 

C  B.  A.  TRENHOLM,  DREA,  JULY  1985. 

C  SOLVES  N  LINEAR  EQUATIONS  IN  N  UNKNOWNS 
C  REDUNDANT  EQUATIONS  PERMITTED 

C  DEFAULT  VALUES  MAY  BE  ASSIGNED  TO  ALL  VARIABLES 
C  ONE  R.H.S.  COLUMN 

DIMENSION  A (MAXN, 1) .X(MAXN) , INDX (MAXN) , JNDX (MAXN) 

C  NOTE  THAT  THE  USER  MASTER  SEGMENT  MUST  CONTAIN  THE  STATEMENT 
C  DIMENSION  A (MAXN, M),X (MAXN), INDX (MAXN), JNDX (MAXN) 

C  WHERE  M  IS  GREATER  THAN  OR  EQUAL  TO  N+l,  AND  WHERE  N  .LE.  MAXN 

C  A.X  *  LAST  COLUMN  OF  MATRIX  A 

C 

C  A(I , 1) .1(1)  +  ...  +  A(I,N).X(N)  -  A(I, N+l) 

C 

C  FOR  I  -  1....N 

C 

C  INDX  IS  A  WORKING  SPACE 

C  JNDX  IS  A  WORKING  SPACE 

C  ISET  -  0  IF  FREE  VARIABLES  ARE  DEFINED  TO  BE  ZERO 
C  ■  1  IF  FREE  VARIABLES  DEFAULT  TO  THE  VALUES  INPUT  IN  VECTOR  X 

C  (FREE  VARIABLES  WILL  OCCUR  IF  THERE  ARE  REDUNDANT 

C  EQUATIONS  OR  IF  THE  SYSTEM  IS  VERY  ILL-CONDITIONED.) 

C  ICK  -  0  IF  A  SOLUTION  X  IS  FOUND 
C  ■  1  IF  EQUATIONS  ARE  INCONSISTENT 

C  (IN  WHICH  CASE  X  IS  RETURNED  AS  A  VECTOR  OF  ZEROS) 

C  IBUG  -  0  TO  SUPPRESS  DEBUG  PRINTOUT 
C  -  1  FOR  DEBUG  PRINTOUT  ON  TERMINAL  UNIT  6 

C  DETERMINE  FLOAT-POINT  ACCURACY,  E.G.  EPS-l.OE-8 
ONE  -  1. 

EPS  ■  1. 

1  EPS  ■  EPS/10. 

EPS1  -  1.  ♦  EPS 
IF (EPS1 . GT . ONE) GO  TO  1 


EPS  -  EPS* 10. 

IF ( IBOG . NE . 0) WRITE (5 . 500) EPS 

500  FORMAT (IX ’EPS  -’E12.5) 

NPLUS1  ■  N+l 

DO  500  I-l.H 

IF (IBUG.NE.O) WRITE(5 ,508)1 , (A(I , J) ,  J-l .NPLUS1) 

608  FORMAT (IX ’ROW '13, ’ :  ’100E12.5) 

600  CONTINUE 

C  SET  INITIAL  SOLUTION  TO  ZERO 
ICK-0 

DO  10  1*1 ,N 

INDX(I)-0 

JNDX(I)*0 

IF(ISET.EQ.0)X(I)*0.0 
10  CONTINUE 

C  NORMALIZE  LINEAR  EQUATIONS  BY  DIVIDING  ALL  COEFFICIENTS  BY  A  CONSTANT 
C  SO  THAT  THE  LARGEST  COEFFICIENT  IS  UNITY. 

EMAX  -  0.0 
DO  16  1*1, N 
DO  15  J-l.N 

EMAX  -  AMAX1 (EMAX, ABS(A(I , J) ) ) 

15  CONTINUE 

DO  16  I-l.N 
DO  16  J*1 .NPLUS1 
A(I, J)*A(I, J)/EMAX 

16  CONTINUE 

C  SOLVE  SYSTEM  OF  LINEAR  EQUATIONS 
DO  20  J*1,N 
ZZ-EPS 

C  SEARCH  COLUMN  J  FOR  LARGEST  ABSOLUTE  COEFFICIENT  >  EPS 
IF(IBUG.NE.0)WRITE(5,501) J 

501  FORMAT (IX* SEARCHING  COLUMN'  13) 

IROW-O 
DO  30  I-l.N 

C  PREVIOUSLY  ASSIGNED  ROWS  ARE  NOT  ELIGIBLE 
IF(INDXU).NE.O)  GO  TO  30 
TEST-ABS(A(I.J)) 

IF(TEST.LE.ZZ)  GO  TO  30 


IROW- I 

30  CONTINUE 

IF ( IBUG . NE . 0) WRITE ( 5 , 5C2) IROW 
502  FORMAT (IX* PIVOT  ROW  -T3) 

IF(IROW.Eq.O)  GO  TO  20 
C  PIVOT  ROW  FOUND 
40  INDX(IROW)-J 
JNDX(J)«IROW 
ZN-A(IROW.J) 

IF ( IBUG . NE . 0) WRITE (5 . 504 ) ZN 
504  FORMATClX’ PIVOT  ELEMENT-1 El 1.5) 

Nl-N+1 

C  NORMALIZE  PIVOT  ROW 
DO  50  K-l.Nl 

50  A ( IROW , K) -A ( IROW , K) /ZN 

C  SWEEP  OUT  OTHER  ROWS ,  FOR  COLUMNS  >  J 
C  (FOR  ALL  COLUMNS.  IF  ISET  NOT  EQUAL  TO  ZERO) 

DO  66  I-l.N 
IF(I.Eq.IROW)  GO  TO  66 
SCALE-A(I.J) 

IF (SCALE. EQ. 0.0)  GO  TO  66 
Jl-J+1 

IF(ISET.NE.O) Jl-1 
DO  60  K-J1.N1 

A ( I . K) -A ( I . K) -SCALE* A ( IROW . K) 

60  CONTINUE 
66  CONTINUE 
20  CONTINUE 

C  CHECK  TO  SEE  IF  UNASSIGNED  ROWS  HAVE  R.H.S.  APPROXIMATELY  ZERO 
DO  80  I-l.N 

IF(INDXU).GT.O)  GO  TO  80 
TEST-ABS ( A ( I . N1 ) ) 

IF (TEST. GT. (EPS* 10.)) GO  TO  99 
80  CONTINUE 

C  MOVE  FREE  VARIABLES  TO  R.H.S.  OF  SYSTEM 
DO  81  J-l.N 

IF(JNDX(J) . GT.O)GO  TO  81 


DO  81  I-l.N 

A(I,N1)»A(I,N1)-X(J)*A(I,J) 

81  CONTINUE 

C  PUT  SOLUTION  VARIABLES  IN  R.H.S.  COLUMN  INTO  VECTOR  X 
DO  70  1*1 ,N 

IF(INDX(I) .EQ.O)  GO  TO  70 
X(INDX(I))-A(I,N1) 

70  CONTINUE 

RETURN 

C  NO  SOLUTION 
99  ICK*1 

IF (IBUG.NE.0)WRITE(5 ,505)1 ,A(I ,N1) 

505  F0RMAT(1X'UNASSIGNED  ROW  #*I4,*  HAS  NONZERO  R.H.S 
C  INCONSISTENT  EQUATIONS 

RETURN 


Appendix  D  Sample  Input  and  Output 

D.l  Terminal  Session  for  Problem  3(a) 


DREA  Tops-20  Command  processor  5. 1(1712) -6 
CRUN  FITBT 
PROGRAM  FITBT 

VERSION  860708  -  DREA  -  B . A . TRENHOLM 
OBSERVED  DATA  FROM  DISK  FILE  LSqiN.DAT: 

IN  ENGLISH  UNITS (0)  OR  METRIC (1)  ?  :  1 
TEMPERATURE  DATA(O)  OR  SOUND  SPEED  DATA(l)  ?  :  0 
READ  HISTORICAL  DATA  FROM  FILE  SPRFO.DAT  (O-NO.  1-YES) :  0 
INPUT  WATER  DEPTH  (  m)  :  200 

OPTIONAL  DEBUG  OUTPUT  TO  TERMINAL  (O-NO.  1-YES)  ?  0 
INITIAL  COORDINATES  OBTAINED  FROM  EXTRAPOLATION 
SUBPROBLEM  1  COMPLETE 

SUBPROBLEM  2  TERMINATES  AS  CASE  I la: EXTRAPOLATED 


LAYER 

* 

1 

CONTAINS 

0 

POINTS 

LAYER 

* 

2 

CONTAINS 

0 

POINTS 

LAYER 

* 

3 

CONTAINS 

1 

POINTS 

LAYER 

* 

4 

CONTAINS 

1 

POINTS 

LAYER 

* 

6 

CONTAINS 

0 

POINTS 

LAYER 

* 

6 

CONTAINS 

0 

POINTS 

LAYER 

* 

7 

CONTAINS 

1 

POINTS 

LAYER 

* 

8 

CONTAINS 

1 

POINTS 

LAYER 

* 

0 

CONTAINS 

1 

POINTS 

LAYER 

* 

10 

CONTAINS 

1 

POINTS 

LAYER 

* 

11 

CONTAINS 

0 

POINTS 

LAYER 

* 

12 

CONTAINS 

0 

POINTS 

2SUBPR0BLEMS  IDENTIFIED 


COORDINATES  OF  THE  FITTED  CURVE  ARE  ON  FILE  FITBT.DAT 
CPU  time  0.77  Elapsed  time  17.06 


D.2  Terminal  Session  for  Problem  3(b) 

CRUN  FITBT 
PROGRAM  FITBT 

VERSION  860708  -  DREA  -  B. A. TRENHOLM 

30 


p: 


r.  -  - 


I 


K- 


.N 


K-v 


K- 

rw-: 


OBSERVED  DATA  FROM  DISK  FILE  LSqiN.DAT: 

IN  ENGLISH  UNITS (0)  OR  METRIC(l)  ?  :  1 
TEMPERATURE  DATA(O)  OR  SOUND  SPEED  DATA(l)  ? 
READ  HISTORICAL  DATA  FROM  FILE  SPRFO.DAT  (O-NO. 

ENGLISH  UNITS  (0)  OR  METRIC  (1)  ?  :  0 
INPUT  WATER  DEPTH  (a)  :  200 
OPTIONAL  DEBUG  OUTPUT  TO  TERMINAL  (O-NO,  1-YES) 
INITIAL  COORDINATES  OBTAINED  FROM  EXTRAPOLATION 


0 

■YES) 


?  0 


SUBPROBLEM  1  COMPLETE 


SUBPROBLEM 

2  TERMINATES 

AS  CASE  lib: CONTROLLED 

9 

LAYER  * 

1 

CONTAINS 

0 

POINTS 

LAYER  # 

2 

CONTAINS 

0 

POINTS 

LAYER  * 

3 

CONTAINS 

1 

POINTS 

m 

LAYER  * 

4 

CONTAINS 

1 

POINTS 

LAYER  # 

5 

CONTAINS 

0 

POINTS 

r  *  7 

LAYER  * 

6 

CONTAINS 

0 

POINTS 

*/  * 

LAYER  # 

7 

CONTAINS 

1 

POINTS 

LAYER  # 

8 

CONTAINS 

1 

POINTS 

i 

LAYER  * 

0 

CONTAINS 

1 

POINTS 

LAYER  # 

10 

CONTAINS 

1 

POINTS 

LAYER  * 

11 

CONTAINS 

0 

POINTS 

»\» 

LAYER  * 

12 

CONTAINS 

0 

POINTS 

■V 

■>: 

2SUBPR0BLEMS  IDENTIFIED 

COORDINATES  OF  THE  FITTED  CURVE  ARE  ON  FILE  FITBT.DAT 
CPU  tlaa  0.89  Elapsed  time  14.78 


i 

D.3 

Input  File  LSQIN.DAT: 

m 

6 

--- 

10.967 

25.3 

9.971 

37.4 

2.938 

80.8 

m 

3.682 

6.653 

93.5 

112.1 

6.024 

137.4 

D.4 

Input  File  SPRFO.DAT 

m 

r  ^  w 

%•  % 

& 

311085 

0 

4600N  6030W  12 

000  45.140  4849.069 

$ 

32 

808  44.780  4847.032 

31 


id 


65.617 

44.240 

4843.685 

98.425 

42.980 

4835.486 

164.042 

39.020 

4808.088 

205.053 

38.210 

4803.349 

246.063 

37.400 

4798.609 

287.074 

37.220 

4798.313 

328.084 

37.040 

4798.017 

410.106 

37.760 

4805.905 

492.126 

38.300 

4811.899 

656.168 

39.380 

4823.798 

D.5  Output  File  FITBT.DAT  for  Problem  3(a) 


12 

16.32845, 

14.20930, 

12.09015, 

9.971000, 

9.971000, 

6.864165, 

3.757331, 

1.991531, 

5.305008, 

6.024000, 

6.024000, 

6.024000, 


0.0000000E+00 

10.00000 

20.00000 

30.00000 

50.00000 

62.60000 

76.00000 

87.50000 

100.0000 

125.0000 

150.0000 

200.0000 


D.6  Output  File  FITBT.DAT  for  Problem  3(b) 

12 

11.94670,  0.0000000E+00 
11.55947,  9.999878 
11.17223,  20.00006 
10.78500,  29.99994 
8.685004,  50.00000 
6.226873,  62.50015 
3.868800,  75.00000 
1.862740,  87.50015 
5.444580,  100.0000 
5.875200,  125.0000 
6.175200,  150.0000 
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